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The Haldane conjecture, when apphed to the Heisenberg 0(3) model with a 6 term in two 
QO . dimensions, states that the correlation length ^ diverges when 9 approaches vr. To verify this 

' conjecture we have numerically simulated the model at imaginary 6 and then analytically continued 

, the results to real 6. We have obtained that the value where the model should become critical is 

CN . 9 = 3.10(5) in agreement with the expectation. 
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I. INTRODUCTION 



It has been shown by Haldane that, depending on the value of the spin cr, the corresponding one dimensional 
antiferromagnetic chains of quantum spins present two kinds of large distance behavior. If a is half-integer, they have 
a power law correlation function. Instead for integer spins they lie in a disordered phase and present an exponentially 
falling correlation function [1,2]. These results were obtained in the limit of large cr. 

The generalization of the above behavior for all values (large or small) of the spin a is called Haldane conjecture. 
[ This conjecture has been widely studied. Actually a partial result had already been proved for a = ^ in [3], while for 
I ■ all half-integer spins it was shown to be correct in [4]. Moreover, the analytic proof for spin a = 1 was given in [5]. 

On the other hand there are indications [2,6,7] that the one dimensional antiferromagnetic chain of quantum spins 
a shares the same large distance physical properties of the two dimensional 0(3) nonlinear sigma model for classical 
spins with a term for 9 = 2Tra. This equivalence would imply that while the ground state of the two dimensional 
"j^ , 0(3) nonlinear sigma model at vanishing 9 must display no long-range order and only short-range spin correlations, 
the model a,t 9 = n should be critical. The first result is well-known, both analytically [8] and numerically [9]. 
i' J However the large distance behavior of the correlation function of the model in the second case is a not so clearly 

■ settled question. 

^ ' The two dimensional 0(3) nonlinear sigma model is a valuable representation of several types of physical problems. 
O Apart from the one dimensional quantum spin chains, in condensed matter physics it may describe the quantum 
Hall effect as well as being useful to understand superconductivity [10]. In particle physics it has in common with 
nonabelian gauge theories some important properties such as instantons, asymptotic freedom (criticality at zero 
temperature), a 9 term, spontaneous generation of mass, etc. 
[ Two recent numerical calculations of the partition function for the 0(3) model in the presence of a term [11,12] 

■ suggest that the theory undergoes a second order phase transition at = tt although the two analyses disagree about 
' the universality class. Indeed the analysis of Bietenholz et al. [11] confirms the critical exponents of the Wess-Zumino- 

Novikov-Witten model at topological coupling fc = 1 as predicted by Zamolodchikov et al. [7], while the numerical 
study of Azcoiti et al. [12] yields a set of continuously varying critical exponents. 

In this work we introduce a direct numerical method to verify the Haldane conjecture for the two dimensional 0(3) 

■ nonlinear sigma model at nonzero 9. The idea is to perform a Monte Carlo simulation to calculate the correlation 
(•~^ length ^ on the lattice as a function of the 9 parameter and to show that it diverges at a precise value of 9, called 

ILJ ^cnd, which, following Haldane, should be 6'end = tt. 
. 5^ ] Due to the (suppossed) divergence of the correlation length at ^ond, a direct simulation would become impracticable 
I as it would require exponentially large lattice sizes. Moreover the Boltzmann weight in the partition function becomes 
^ ■ complex for real 9 and consequently it loses its probability meaning, thus precluding the importance sampling of Monte 
Carlo methods. We overcome these two difficulties by simulating the theory at imaginary 9 (where ^ turns out to be 
small enough to allow the use of moderate lattice sizes) and analytically continuing the results to the real 9 values. To 
this end we introduce a new fast cluster algorithm that works for imaginary nonzero theta. This work is an extended 
version of the paper appeared in [13]. 

In the next section we shall discuss the formulation of the model on the lattice and the corresponding lattice 
definition of the 9 term and its meaning. The method to calculate the 9 term is introduced in Section HI. The 
new cluster algorithm expressly devised for the present work shall be described in Section IV. The results and 
corresponding plots are displayed in Section V. We end the paper with some conclusive comments in Section VI. 
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II. LATTICE IMPLEMENTATION AND NOTATION 

The Boltzmann weight of the partition function in our simulations was exp {—S) with 

S = A-ieQ, 

where the sum is extended over all lattice sites x and directions /i = 1,2. The factor (3 is the inverse temperature (in 
units of the spin coupling), Q is the total topological charge or winding number of the configuration (see later) and 
(l){x) is a 3-component unit vector that represents the dynamical variable, a classical spin, at the site x. We have 
used a square lattice of lateral size L with periodic boundary conditions. 

In the limit where the lattice spacing a vanishes (keeping L ■ a fixed) , the above expression for A becomes the action 
of the classical field theory defined on a continuum two dimensional plane 




together with the condition cf){x)'^ = 1. 

In Fig. 1 we show a stereographic projection that defines an instanton configuration on the 0(3) nonlinear sigma 
model. The two dimensional plane is shown where the configuration of spins lies. There is a unit sphere resting on 
the origin of the plane. One can draw straight lines that join the north pole N of the sphere with an arbitrary point 
F on the plane. One such a line pierces the sphere surface at point P. The unit vector that begins at the center C 
of the sphere and points to P is the value of the spin vector 4>{F) to be assigned at the point F of the plane. This 
construction defines a configuration called instanton and its analytical expression is the following (the lattice has been 
replaced by a continuum by sending a ^ as above, furthermore polar coordinates r and (f are used to locate the 
position of the spin variable) 

7*. , /4rcos(^ irsiiKf — 4\ 

In particular, notice that all spins at infinity are identified with the same value (j){r = oo,ip) = (0,0, 1). Then spins 
slowly rotate while approaching r = until becoming (j){r = 0, t^) = (0, 0,-1) at the origin. 

The spin value in Eq.(3) is a solution of the classical field equations associated to the action Eq.(2). 



N 




FIG. 1. A stereographic projection defines an instanton on the two dimensional 0(3) nonlinear sigma model. 

The main feature of instanton configurations is that the set of all spin vectors describes a complete winding of the 
sphere, as it is obvious in the example shown above. This winding can be calculated by the integral [14] (valid in the 
continuum two dimensional plane) 

Q = y d^xQ{x) , 

Q{x) = ^e^'-'ebcdct>\x)d^ct>%x)d,cl>''{x) , (4) 
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where Q is called topological charge or winding number and Q{x) is the topological charge density. Spatial indices 
H, V and 0(3) vector indices 6, c, d are summed up. For the instanton of Eq.(3) it yields —1. There arc however many 
more instantonic configurations, besides that shown in Fig. 1 and in general Q takes any positive, negative or null 
integer value, depending on how many times and in what direction the whole set of spins covers the unit sphere. 

The second main property is that instanton configurations carry a finite amount of energy (the operator A in Eq.(2) 
takes on a finite value) even when the lattice size L diverges. 

In the Monte Carlo simulations we have used two different definitions of Q{x). The first one [15] 

(y^\x + V)-<\><'{x-v)y (5) 

is a symmetrical discretization of the expression for Q{x) in Eq.(4) and is usually called "naive" definition. The 
corresponding winding number is Q^^^ = Q^'^\x). 

The second lattice expression that we used in our simulations is defined on triangles (not on single sites). Every 
plaquette of a square lattice can be cut through a diagonal into two triangles. If we call ^i, (^^ and 4>z the fields at the 
sites of the three vertices (numbered counterclockwise) of one of these triangles, then the fraction of spherical angle 
subtended by these fields is Q^^^A) and it satisfies [16] 

exp (27rzQ(2) (A)) = ^ (l + • ^2 + • 03 + ^3 • <?i 

• ((?2 X (^Ts) ) , (6) 

where (? = 2(l + (/)i •(/)2)(l + (/>2 •03)(l + 03 -^i) and (5*^^)(A) G [— ^, +5]- The above conditions uniquely determine the 
portion of spherical angle subtended by 0i, ^2 and ^3 and the sum of (5'^^(A) over all triangles yields the so-called 
"geometrical" topological charge Q^^^. 
The two definitions Q^^'^^ belong to the same universality class. In particular they both satisfy the limit 

Q^^^^\xf-^a'Q{x) , (7) 

where a is the lattice spacing. 



III. EVALUATION OF Q 

In general, a definition of Q on the lattice does not necessarily lead to integer values on a single configuration. 

To recover integer results for Q on ensembles of configurations that contain the same topological charge, we must 
renormalize this operator. The lattice and the continuum topological charges are related by [17] 

Q(i.2)^^a.2)g^ (8) 

Zq ' being the corresponding renormalization constant. The origin of this constant can be traced back to the 

presence of statistical fluctuations in an otherwise smooth instantonic configuration. In general operators that reveal 

the topological charge Q of a configuration give wrong answers due to the disturbance caused by the presence of 

fluctuations. 

(12). 

The Zq ' function depends only on the temperature /3 and is chosen in such a way not to depend on since the 

(1 2) 

introduction of this term does not modify the structure of fluctuations in the model. Moreover it satisfies < Zq' ' < 1 
for all values of the lattice spacing [18,15]. 

Zq''^^ can be calculated either in perturbation theory [17,18,15] or by a nonperturbative numerical method [19-21]. 

We have used the latter. In a nutshell it works in the following way: one measures Q*^^'^^ on an instantonic configuration 
(topological charge Q = +1) after heating it at a temperature f]. From Eq.(8) this measurement yields Zq' 



First of all an instanton with topological charge +1 is put by hand on the lattice. We used the solution [22] 

+ i4>^ L/2 - i{x2 - L/2) 

~ AeW4 ' 



(9) 
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where (jf' is the &-component of the field at site x = {xi,X2), L is the size of the lattice and A is the size of the 
instanton. This solution is the one shown in Eq.(3) after centering it in the middle of the lattice (L/2, L/2), dilating it 
to the size A (the size in Eq.(3) is A = 2) and performing appropriate rotations both in a; space and 0(3) space [23]. In 
order to work with a rather stable instanton (recall that single instantons on the lattice are only metastable solutions 
of the equations of motion on a torus) it is convenient to choose A < 0.15 L and A > 8 [21]. We took A = 16 on a 
L = 120 lattice. 

Then 100 updating steps are applied (we used the Heat-Bath algorithm [24] on the conventional 0(3) nonlinear 
sigma model without a term since the renormalization constant to be used in Eq.(8) does not depend on 0). After 
every Heat Bath step the value of Q^^^^) jg measured and, in order to monitor the instantonic contents and check that 
it is not varied after the updating step, Q^^'^-* is measured again after 6 relaxation hits applied on a separate copy of 
the running configuration. The complete history of 100 Heat-Bath updating steps and related measurements of Q^^'^^ 
is called a trajectory. In the calculation of Zq^ we used 4 • IC* trajectories at /? = 1.5 and 1.6 and 10^ trajectories for 
/? = 1.7 and 1.75. The average of Q^^'^) on all trajectories, as long as their topological charge remained equal to +1, 
yielded Z^q'^\ 

The 0(3) nonlinear sigma model develops an infrared divergence in its instanton size distribution [25,20]. This 

divergence facilitates the copious creation of new instantonic objects at every updating step, thus modifying the total 
topological charge of the configuration. For this reason the relaxation test is extremely important. 

As a relaxation method we used the so-called cooling [26]. It sweeps through the whole lattice and modifies one 
by one every single spin variable in order to locally minimize the energy of the relax;ed configuration. Actually many 
variants of the cooling method exist in the literature but in practice all of them act analogously [27]. 

The above nonperturbative method is summarized by the expression 

7(1.2) _ /l-instanton^'^Q^'''^ ^xp (-^) 

— n zz^ ; — n — ' \^^) 



where Vcj) stands for the measure 



Xl— instanton '^'^ ( 



l[is{${xr-i) nd/(^) (11) 



and /i_instanton means that the integration is extended over all configurations (fluctuations) that preserve the back- 

groimd of one instanton. Since the geometrical charge Q^^) jg j^i till the backgroimd configuration is one instanton, 

(2) 

the expression (10) yields Zq = 1 for all /3 [28]. This result derives from the fact that fluctuations, viewed as local 
large (positive or negative) values of the spherical angle in some spherical triangles, cancel out when summing up all 

individual contributions Q(2)(A). 

The determination of Zq ^ is not so trivial and an example of such an evaluation is shown in Fig. 2. Measures of 
Q^i) on configurations that have topological charge +1 attain to a plateau (in general after a few Heat-Bath steps) 
and stay on it for the rest of the updating steps. The height of this plateau is the value of Zq\ In Table 1 the results 
for Zq^ at the values of /3 used in the present work are given. 



Table 1. Zq and ^end for the topological charge Q'^^K 



/3 


(^^cnd) 




X^/d.o.f. 


^'ond 


1.5 


111(5) 


0.285(9) 


0.90 


3.00(12) 


1.6 


94(5) 


0.325 (G) 


0.45 


3.15(10) 


1.7 


67(3) 


0.380(6) 


1.04 


3.11(9) 


1.75 


56(3) 


0.412(5) 


0.68 


3.08(9) 
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FIG. 2. Data for (Q'^') start at +1 at the 0-th Heat-Bath step and then they go down until reaching a plateau. The height 
of the horizontal line and grey band are the value and error respectively of ZQ\f3 = 1.5). 

In Fig. 3 an histogram of the distribution of topological charge Q*-^' is shown. It has been produced from the data 
of Q'"^'' obtained during the calculation of Zq ^ at /? = 1.5 and contains 70 bins within the interval Q^^^ G [—4, +4]. 

For each trajectory in the calculation of Zq^ we obtained the average of Q^^^ over the steps that come after the onset 
of the plateau (from Fig. 2 this happens at about the 15th step for /3 = 1.5) and for which the cooling test gave a 
background topological charge +1. The histogram of Fig. 3 displays the distribution of these averages. Each of the 
above averages turns out to be an uncorrelated estimate of Zq K The thick vertical hne is the value of Zq \ 0.285(9). 

Actually the error in the evaluation of Zq^ was determined by using this kind of plot. Indeed, it was extracted by 
usual gaussian analysis on the histogram. 

Recall that the number of trajectories for (3 = 1.5 was 40000. However the area of the histogram in Fig. 3 is much 
less than 40000. This is due to the fact that in many trajectories the background charge is no longer +1 already at 
the beginning of the plateau. In such cases, the whole trajectory was discarded. The area under the histogram in 
Fig. 3 is about 6700. The small ratio 6700/40000 gives an idea of the frequent creation of new instantonic objects 
that modify the background topological charge. 
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FIG. 3. Histogram displaying the distribution of naive topological charge at /3 = 1.5 starting from an initial configuration 
with background topological charge +1. Observe the shift of the gaussian-like distribution towards positive values of Q^^^ (a 
vertical thin line indicates the zero value and a thicker line the result of Zq^). 

The; renormalizatioii of the topological charge brings about a relevant consequence for our study: the parameter 
that appears in the expression of the Hamiltonian used in the computer program during the simulations in general is 
not equal to the true physical 6 parameter. Henceforth we shall call 9 the parameter that appears in the simulation 
program and the relation among the two parameters is ^ = 6Zq ' . The value of B where the correlation length 
diverges will be called ^end- Since Z^q = 1, it is clear that this distinction among theta parameters is irrelevant for 
the geometrical charge Q'^^ . 
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Wc have simulated the model at several temperatures f3 and parameters 9 for the two topological charge operators 
Q*-^-* and Q*-^-* in order to show that our results are independent of the operator chosen for the simulation. Moreover, 
for the case of the naive topological charge Q*-^-*, we have introduced a very fast cluster algorithm. Instead, for the 
geometrical charge Q^"^^ a rather slow Metropolis-like algorithm has been used. Actually the cluster algorithm for the 
naive charge was so efficient that it took much less computer time to investigate this charge than the geometrical one 
even though the naive charge required the extra calculation of Zq^ by a separate out-of-equilibrium simulation for 
each temperature /3 as described above. 



IV. CLUSTER ALGORITHM FOR IMAGINARY 9 

Although the use of the topological charge density Q^^^ requires the knowledge of a renormalization constant, it 
brings about the advantage that the Hamiltonian in (1) can be simulated on the lattice by use of a fast cluster 
algorithm when 9 is imaginary. 

Let us briefly describe the main characteristics of the new cluster algorithm expressly devised for the present work. 
The first part of an updating step with the usual Wolff algorithm [29] for the standard 0(3) sigma model without a 9 
term consists in choosing a random unit vector r in such a way that every dynamical field can be split in a component 

parallel to r and the rest, ${x) = (^${x) ■ r + $±{x), where $±{x) denotes the part of ${x) orthogonal to r. Then 
the signs of (^${x) ■ for all x are updated a la Swendsen-Wang as in the Ising model [30] . 

By introducing the above separation for <f>{x) in the expression (5) and recalling elementary properties of determi- 
nants in three dimensional vector spaces, we can rewrite it as 



QW(,t) = 
1 

167r 



I (${x) ■ f) (di,2 + rf-1,-2 + ^2,-1 + C^-2,l) 

+ ((^(a;+i) -r) (do,-2 - ^0,2) 
+ ((?(x-T)-r) (do,2"do,-2) 
+ ((^(x + 2) -r) (do,i - c«o,-i) 

+ (<?(x-2).f) (do,-i-tio,i)} , (12) 

where a; ± 1 means the site at the position one step forward (backward) in the direction "1" starting from site x and 
the notation dij stands for the 3x3 determinant 

(f)\x+7) (t)^{x+l) (t>^{x+l) . (13) 

(f)^{x+T} (t>^{x+3) (t>^{x+3)/ 

In this fashion the theory at each updating step looks like an Ising model in the bosom of an external local magnetic 

field h{x) because the expression in Eq.(12) is linear in • . The value of this field varies at each updating and 

accordingly it must be recalculated after every step. Recall that all Monte Carlo simulations have been performed 
with an imaginary parameter 9 = E R). By gathering all contributions of the type shown in Eq.(12) that 

contain (j){x) ■ at site x one can readily derive the effective magnetic field at this site, 

^ - / 

h{x) = - Y^l«^(a;) ■r\[di^2 + rf-1,-2 + ^2-1 + d-2,1 

+ C^-1,-1-2 + rf-1+2,-1 + rfl,l+2 + di-2,1 

+ ^2,2-1 + ^2+1,2 + cJ-2,-2+1 + d-2-1,-2 ) • (14) 



di+k,j (and analogous terms in (14)) arc the straightforward generalization of the above definition (13) when the site is 
obtained by shifting two steps from the original position a;, the first in the direction i and the second in the direction 
k. 
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Hence the last step in the updating consists in applying to the above expressions an algorithm valid for the Ising 
model in presence of a magnetic field. In the literature there are two such algorithms, the Lauwers Rittenberg [31] 
and the Wang [32,33] methods. After testing their perfomances and comparing the corresponding decorrelation times 
with the usual Metropolis [34] , Heat-Bath [24] and overHeat-Bath [35] , we decided on the Wang algorithm. It consists 
in placing the magnetic field on an extra, fictitious site (called ghost site or ghost spin) that couples to every Ising 
spin through the value of h{x). Using this coupling on the same footing as all other terms in the action, the Fortuin- 
Kasteleyn clusters [36] are created by using the Hoshen-Kopclman algorithm [37] and then updated with the usual ^ 
probability. The only distinctive feature of the presence of a ghost spin is that the cluster that contains it does not 
flip. 




(1) 



Cluster 2 
Metropolis Q' 



(2) 



FIG. 4. Autocorrelation functions for the Metropolis (Q''^') 
time" T which has a discrete ticking at each updating step. 



80 

1(1) 



100 



and cluster {Q^ ') algorithms as a function of the "updating 



Following the proof given in [29] , it can be seen that our algorithm also satisfies the detailed balance property. 
To generate the initial random unit vector r the method proposed in Ref. [38] was used. 
In Fig. 4 we show the autocorrelation functions 



C{r) ^ ^^"^ 



/ T7\ 2 



(15) 



calculated from the measures of the energy operator E = (p{x) ■ (l){x + 'fl) (not summed over /x), for the two algorithms: 
Metropolis when the Q*^^^ operator is used and the above described cluster algorithm for Q^^h Er indicates the r-th 
measurement. In both cases f3 = 1.5 and L = 120. The theta parameter was ■& = 10 for Q''^'' and ■& = 2.85 for 
Q'^^ (note that this choice was dictated by the condition 6 = Zq\i3 = 1.5)6). The plot clearly exhibits the major 
efiiciency of the cluster algorithm. 



V. RESULTS 



As the ground state of the model is a triplet [39], we studied the correlation functions of operators having one 0(3) 
index. We measured the correlation of the two operators 

^i{x) = ${x), ^2{x) = ${x)x${x + l). (16) 

Firstly we calculated the corresponding wall operators by averaging over the xi coordinate 

M)i(X2) = -^^^i(x), y^.ix,) ^ 1^^2(X) . (17) 

Xl Xl 

In order to extract the correct correlation length and to clean its signal from any mixture with higher eigenvalues of 
the Transfer Matrix, we used the variational method of Ref. [40], where ^ is obtained from the exponential decay of 
the largest eigenvalue of the correlation matrix 
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(Wi(x2) • Wjm - {Wi) ■ (W,) . (18) 

No improved estimators [41] were used since the operators in (17) contain too many fields <!) and this fact leads to 
intractable sums over clusters [42]. 




FIG. 5. Determination of ^ at /3 = 1.5 for four values of Im{6) = •& for the case of the Q^^^ operator. Lines are drawn to 
guide the eye. 



In Fig. 5 the method that we followed to extract ^ is illustrated with four examples. The exponent of the fall-off of 
the largest eigenvalue amax of the correlation matrix, Eq.(18), is plotted against the distance X2- For each X2 it was 
extracted by comparing the expression 



Ctr, 



Q!max(a;2 - 1) 



(19) 



with the theoretical behavior {L is the lattice size) 



cosh ( (X2 - L/2) /e 
cosh ( (a;2 - 1 - L/2) 



(20) 



and its error was determined by jackknife. An approximate plateau is clearly identified at moderate distances. The 
definite values of ^ and its error were chosen self consistently at X2 — 2^. 

For each /? we calculated ^ from Monte Carlo simulations at several imaginary values of 6 for both Q^^' and Q^^' . The 
set of results for a given (3 were then analytically continued from imaginary 9 to the real axis by an usual procedure 
of numerical extrapolation. In the extrapolation we avoided using a trial function dictated by some theoretical 

(— 2\ 2/'^ 
C2 — 9 j which is, up to logarithmic corrections, the Renormalization Group 

prediction [43], because such an analytic form implicitly assumes the vanishing of 1/^ at a precise value of 9 (not 
to say that it is supposed to be accurate only in a close neighborhood of its zero, 9 = ^foi). Instead we made the 

2 

extrapolations by using polynomials in 9 and ratios of such polynomials. These functional forms are indeed both 
simple and very general and they leave room for any possible behavior in 9. 



8 



FIG. 6. Four extrapolations for the same set of Monte Carlo data of 1/^. The scale in the axes is the same for the four 
windows. Data were extracted from simulations at = 1.5 and using the Q^^-* operator for the topological charge. The 
horizontal black bar is placed at the value of 9 where the Haldane conjecture predicts a critical behavior. Each continuous line 
is the result of the extrapolation and the dashed lines enclose the boundary of its error. 



A. Results for Q'-^^ 

2 ■ 10^ decorrelated propagators were measured for all values of 9 at each /3. They were obtained after separating 
consecutive configurations by a combination of one Heat-Bath, two overHeat-Bath and one cluster updatings. The 
values obtained for 1/^ at /3 = 1.5 with their error bars are the squares in Fig. 6. In this figure four different 
extrapolations are shown (the extrapolation functional forms are displayed). The Haldane conjecture predicts that 

1/^ vanishes at ^7r/Zg^(/3 = 1.5)^ . This value is indicated by the small horizontal shadowed bar on the 0^ axis (its 

horizontal width arises from the error in the evaluation of Zq \i3 — 1.5), see Table 1). Prom Fig. 6 it seems clear that 
all analytic continuations arc in fair agreement among themselves and with the Haldane conjecture. We emphasize 
that no prejudices about the possible zeroes were included in the extrapolating functions. 
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FIG. 7. Behavior of 1/^ as a function of 9 . Circles {(3 = 1.5), up triangles (/3 — 1.6), squares (/3 = 1.7) and down triangles 
(/3 = 1.75) are the data from the simulation at imaginary 9 {9 < 0) by using the Q^^^ lattice topological charge. Meaning of 
continuous and dashed lines as in Fig. 6. 

Four values of /? were studied in the case of Q^^^ : (3 = 1.5, 1.6, 1.7 and 1.75. The respective sets of Monte Carlo data 

for 1/^ are shown in Fig. 7 as circles, up triangles, squares and down triangles. The lattice sizes were 120, 180, 340 

2 2 

and 470 respectively. The extrapolations in this figure were done by using the functional form (ci + C2 ^ )/(l + C3 ^ ) 
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since it was the one that produced the value of the statistical test closer to 1 for all four (3's (the values for x^/d.o.f. 
are listed in Table 1). 

The results of the analytic continuations are given in Table 1. The physical value of 6 where the theory becomes 

critical is given by ^cnd = ^^ond-^Q''- The numbers in the last column are in good agreement with the prediction 
that criticality is achieved when 9 equals tt. Other functional forms used for the extrapolations led to similar results 
although in some cases the test was far from unity and hence the related extrapolation seemed statistically unlikely 

_2 4 _2 

(sec for instance the (ci + + cr^B + c^O ) case in Fig. 6). 

The lattice sizes were chosen large enough to meet at 6* = the condition L/^ > 10 (to be specific, these ratios 
were 10.8, 9.5, 9.8 and 9.9 for /3 = 1.5, 1.6, 1.7 and 1.75 respectively). Once this inequality holds at ^ = 0, it is 
amply realized at the values of 6 where the simulations were performed as inferred from Fig. 7. This fact warrants 
the absence of significant finite size effects. 

On the other hand, from the Monte Carlo values of ^ shown in Fig. 7 one can see that we simulated the model at 
correlation lengths that altogether satisfy ^ > 5 which is safely far from the strong coupling region (in this region the 
universality property loses its meaning and results may depend on the choice of operators used in the action). 

B. Results for Q*^) 

In this case the visual Metropolis algorithm was used for updating. Observe that when the spin variable 4'{x) is 
updated, the 9 term in Eq.(l) will contribute to the variation of the Hamiltonian only if the topological charge gets 
modified within the only six triangles that surround the site x. Such modifications occur only if some instantonic 
object rises or disappears in the area delimited by these six triangles. Such an event barely occurs on a such a small 
area and as a consequence prolonged decorrelations must separate consecutive measurements of any operator (see 
Fig. 4). 

10^ independent propagators, separated by 100 decorrelation updatings, were measured for each value of 9 (recall 
that for Q^'^^ we have 9 = 9). We report data for only two values of (3. Notice that the total statistics and the number 
of /? values studied is evidently smaller here than in the previous subsection. As explained above, this is due to the 
use of a much less efficient updating algorithm. 

Data are displayed in Fig. 8 as squares and triangles for ^ = 1.5 and 1.55 respectively. The corresponding lattice 
sizes were L = 110 and 150. The value of 9 where Haldane predicted the vanishing of 1/^ is indicated with an arrow, 
9^ = TT^. Data near 9^ = are very noisy. This is due to the relatively poor statistics obtained in the simulations 
with the operator Q^'^\ 

The numerical results are given in Table 2. Comments similar to the Q^^-* case apply to the extrapolations shown 
in the figure. From Fig. 8 the correlation lengths satisfy ^ > 3 which is still away from the strong coupling regime. 
Again the results are in fair agreement with the conjecture. 



Table 2. ^gnd for the operator Q^'^\ 



/3 


(^^cnd) 


7(2) 


xVd.o.f. 


^^cnd 


1.5 


10.4(1.0) 


1.0 


1.72 


3.22(16) 


1.55 


9.7(1.0) 


1.0 


0.73 


3.11(16) 



By averaging all results for both topological charge operators and assuming gaussian errors we obtain that the 
model should become critical at ^end = 3.10(5). This is the chief result of our work. 
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-12 -6 6 12 

FIG. 8. Behavior of l/£, as a function of 6^ for the geometrical charge Q^^K Squares (/3 = 1.5) and up triangles (/3 = 1.55). 
In this case 6 = 6 and the position of ^ = tt is marked. Meaning of continuous and dashed lines as in Fig. 6. 

VI. CONCLUSIONS 

We have simulated the 0(3) nonhnear sigma model in two dimensions with an imaginary term at several values 

of the inverse temperature f3. The correlation length was measured and extrapolated towards real 9. In all cases 
the extrapolation showed a divergence at a value of 9 compatible with the Haldane conjecture = n. Our result 
is ^ = 3.10(5) which agrees within errors with the conjecture. This value seems very robust as it is independent of 

the topological charge density operator chosen for the simulation. In particular, an operator Q^^"^ that requires a 
nontrivial renormalization constant leads to the same conclusion than another operator (the geometrical charge Q^^-*) 
that does not renormalize. 

A direct numerical study of the model at 9 values which are both real and close to tt is unfeasible. First of all at real 
9 the Hamiltonian becomes complex, thus preventing the importance sampling in Monte Carlo methods. Furthermore 
exponentially large lattice sizes would be required to avoid the severe finite size effects that would supervene in 
simulations performed close to the critical point 9 = n. 

A new fast cluster algorithm was purposely introduced to simulate the theory with an imaginary 9 term. It works 
for the operator Q^^h 

A salient outcome of our work is the good performance of the analytic continuation from imaginary to real 9. No 
theoretical prejudices were assumed in the functional form used in the extrapolation and different functions led to 
comparable results (see Fig. 6). 



10 

8 

Im(Q) ^ 

4 
2 



12 3 

P 

FIG. 9. Possible phase diagram in the Im{6)-l3 plane, {Im{d) = t9). The phase transition lines (continuous line is first 
order and dashed line is at most second order) are borrowed from Ref . [44] . Symbols are situated at the positions where single 
simulations were done. Filled (white) symbols mean simulations with Q'^^ (Q'^^) and the vertical sequences of data correspond, 
from loft to right, to (3 — 1.5, 1.55, 1.6, 1.7 and 1.75. In this diagram extrapolations appear as going downwards starting from 
the Monte Carlo data points, hence they never touch the phase transition lines. 
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A key ingredient for the successful extrapolation was to have got data from simulations within a wide range of 
(imaginary) values of 6 for all /3, = -iO e [0, 10] when Q^^' was used and d ^ d = -iO € [0,3.5] for Q*^^ the 
difference of intervals being due to the effect of the nontrivial renormalization that must be applied to the former). We 
also noticed that at very large values of imaginary 9 the extraction of the correlation length becomes more imprecise 
since the kind of plateaux shown in Fig. 5 shrink making it difficult to decide the correct result for ^ and its error 
bar. This problem appears approximately at ?9 > 15 for the Q^^-* charge density. Then, the method consists in taking 
data in a wide range of values of imaginary 9 [ov 9) while staying not too far from the real axis in order to make the 
extrapolation sensitive to the real physics of the problem and in order to have wide enough plateaux which allow to 
easily extract the correlation length ^. 

Moreover our extrapolations did not cross any nonanalytical region in the phase diagram of the model. Indeed in 
Fig. 9 we show a schematic reproduction of the phase diagram that appears in Ref. [44]. The first order transitions 
lie on the continuous straight line while the dashed line indicates less singular transitions, at most second order. Our 
simulations were performed at the positions indicated by the symbols (same symbols than in Figs. 7 and 8; notice 
that this plot refers to the values of 9, not of 0, for both topological charge density definitions) and the corresponding 
analytic continuations proceed downwards, thus very far from any possible line of singular points. 

For 9 > IT the model should acquire again a finite correlation length since the Hamiltonian is a periodic function of 6. 
However an analytical continuation could hardly display such a behavior beyond 9 = it because of the nonanalyticity 
at that value of 9. 
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